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Abstract 

At high gluon or string densities, gluons' saturation or the strong interaction 
among strings, either forming colour ropes or giving rise to string's percola- 
tion, induces a strong suppression in the particle multiplicities produced at 
high energy. This suppression implies important modifications on cosmic ray 
shower development. In particular, it is shown that it affects the depth of 
maximum, the elongation rate, and the behaviour of the number of muons at 
energies around 10 17 -10 18 eV. The existing cosmic ray data point out in the 
same direction. 

One of the most crucial astrophysical issues of the highest energy cosmic rays (above 
~ 10 17 eV) is that of their composition. This problem is linked to the identification of 
the origin and possible sources of these cosmic rays. Current theoretical models expect a 
transition from galactic to extragalactic or galactic halo sources near the region of the ankle 
which leads to the usual expectation of the changing of composition from heavy to light 
elements. 

Experimentally, measuring the composition at these energies is a challenging task. The 
very low fluxes involved imply that one has to rely on indirect measurements which depend 
on simulations of the development of cosmic ray cascades in the atmosphere. These, in 
turn, are model dependent and, specifically, depend on extrapolations of hadronic models 
to energies and kinematical regions never measured in the laboratory. There is, therefore, 



some degree of uncertainty in the shower development and one may ask what is the effect of 
this uncertainty in the reconstruction of shower parameters, mainly total energy and mass 
composition. 

To avoid this problem, experimental groups have concentrated on observables which 
are expected to be more or less independent of the hadronic model used, or which have 
its dependence under theoretical control. These parameters include the maximum of the 
cascade, A max , the slope parameter, (3 = dlog(p fi (6Q0)) / d\og{E) , where ^(600) is the muon 
density at 600 m from the core, and the elongation rate, D w = dX miLX / d\og 10 (E) . Other 
parameters have been less frequently used, see Ref. J7| for a general review. 

Several experimental groups have measured the cosmic ray spectrum and mass compo- 
sition in the ankle region and beyond using the above mentioned parameters [0-0 see also 
0. The results on mass composition are inconclusive. Fly's Eye observe a change on the 
slope of X max in the region around 3 x 10 17 eV, which is interpreted as a change on the com- 
position from heavy (iron) dominated to light (proton) dominated. However, AGASA |2j 
measures a muon component and (3 parameter consistent with iron on that region. Although 
some part of the discrepancy between AGASA and Fly's Eye may be due to the use of dif- 
ferent hadronic models ||, as pointed out by Nagano et al. |J the important issue is that 
AGASA sees no significant change on the muon component of showers all along the ankle 
region and beyond, from 10 16 ' 5 eV to 10 19 ' 5 eV, and thus no strong change on composition 
is inferred. HIRES and MIA collaborations 0] have jointly measured both the X max and 
the /3 parameter. They observe a strong change of X msx with energy, which implies a large 
elongation rate, D w = 95 gr/cm 2 and on the other hand they see no change on the slope of 
the muonic component, j3 = 0.73, measured which is broadly compatible with the AGASA 
observations. HIRES and MIA however have measured these parameters in a narrow range 
of energy, from 10 17 to 10 18 eV and with low statistics, only during a limited exposure. 

In this paper we will show that under rather general conditions a change on the hadronic 
interactions at the energies of interest is expected, which may have important consequences 
for the interpretation of cosmic ray data. Whether this change is enough or not to produce 
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the observed changes on the cosmic ray data we can not tell at present. On the other 
hand we can state the necessary conditions for this change to explain the observed data: 
i) There should be an abrupt change on the hadronic interactions at the observed energy 
-Slab ~ 5 x 10 17 eV for Fe-Air collisions. This corresponds to a CM energy of ~ 4200 AGeV 
and a density of gluons of ~ 9 fm -3 ii) At this energy the slope of the growth of the 
multiplicity with energy should vary from ~ 0.24 to a essentially flat < 0.09. If i) and ii) are 
verified then there is no additional need for a change on composition to explain the data. 

It is important to point out that, although the change on the multiplicity may or may 
not be enough to produce the observed results, some effect should always be present and 
should be taken into account in any realistic simulation of cosmic ray showers. Currently 
no Monte Carlo code for cosmic ray showers has yet been implemented with these effects^ 

In the last years, a wealth of data coming from HERA and the heavy ion SPS experiments 
have risen questions about the behaviour of the hadronic interactions at very high energy. 
We may consider perturbative, gluons, or non perturbative, strings, as the fundamental 
variables of our description, At high gluon density, the saturation of gluons |9| and/or a 
strong jet shadowing [[TtJ are expected. In the case of high string density we expect the 
fusion of strings [TTJ or colour rope |T2] |, and probably, above a critical string density the 



percolation of strings |L3[ and the formation of quark gluon plasma at the nuclear scale are 
expected. 

One general feature of all these hadronic phenomena is the strong suppression of particle 
multiplicity compared to the multiplicity expected in their absence. Namely, for central 
Pb-Pb collisions the charged particle multiplicity expectations in the central rapidity region 
changes between 1500 (7500) for the relativistic heavy ion collider, RHIC, (the large hadron 
collider, LHC) for models that do not include these effects to 900 (3000) when they are 



*Sibyll version 2.0 § incorporates some shadowing effect. However this was done for pp collisions 
only and does not affect to our reasoning below. 
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included [14]. As a framework we will use the quark gluon string model (QGS) ||15|| , a 



modified version of the Dual Parton Model |16[] . The model is based on the large N expansion 
of QCD but it is largely phenomenological and describes most of the soft hadronic physics 
rahter well. Inclusion of hard, perturbative, physics has been done in various ways. In the 
quark-gluon string model, multiparticle production is related to the interchange of multiple 
strings which break and subsequently hadronize. 

In this model one can most easily understand the expected changes on the behavior of 
hadronic collisions at high energies. It is more convenient to work in the plane transverse 
to the collision. In this plane, strings are seen as small circles of fixed radius, r. As the 
energy increases, the number of strings interchanged increases and the total area occupied by 
strings increase. At high energy, strings start to overlap and fuse together. For high enough 
string density, n c , strings may percolate in a second order phase transition, i.e. continuous 
paths of strings are formed in the collision area. Since the number of independent strings is 
reduced after the fusion one expects a depletion on the number of particles produced, i.e. a 
reduction on the multiplicity. 

In the QGS the multiplicity grows with energy as n(s) ~ s A , where A is related to the 
intercept of the soft pomeron |]T3|[|. In the case of percolation, the reduction of multiplicity 
is given by fli~7,14| 



n'(a)= nW^)), (1) 

where 

Hv) = l ~ Z f 1 , (2) 

and the parameter rj is the fraction of the total area occupied by strings 

nr 2 N s , . 



^Here we consider minimum bias events, which are relevant for cosmic ray experiments. The 
parameter A may depend on the centrality of the collision 
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Here N s is the number of strings produced in the collision, r is the string's transverse size, 
and R is the total collision area. N s grows with energy as iV s ~ s A ' , where A' is the intercept 



of the soft pomeron ]nj. Therefore at large rj the total multiplicity grows with energy as 

n'(s) ~ s A ~ A '/ 2 . (4) 

This reduction of multiplicity is not exclusive of the percolation of strings. For instance, 
in perturbative QCD a reduction in the number of jets produced as the energy increases 
is also expected [fLOR . At high energy the number of jets produced grows with energy as 
n(s,Pt) ~ (s / (4p 2 )) AH , where Ah is the intercept of the hard pomeron and pt is the trans- 
verse momentum of the jet. A (mini)jet occupies a transverse area of order ft/pf, since the 
number of jets increases rapidly with energy, saturation occurs when the area occupied by 
the jets equals the total transverse area [|I~0| , |14]1 : 



which implies 

n(s,p 2 t ) ~ . (6) 

Both Eqs.(f|J|) have been checked directly in Monte Carlo simulation [p] , ^0[| . Surpris- 
ingly, for nucleus-nucleus collisions the reduction in the power of multiplicity growth with 
energy is of the same order both for the case of string fusion and of shadowing. The power 
changes from ~ 0.24 to ~ 0.19. In general parton saturation, shadowing, string fusion, or 
percolation will produce the effect of reduction of the multiplicity although we expect the 
degree of this reduction to be model dependent. 

For cosmic ray showers, the rate of change of the multiplicity with energy is directly 
related to the change of the elongation rate. This has been known for a long time as the 



elongation rate theorem ||18|| . The elongation rate theorem can be deduced easily, as it follows 
from a scaling argument. Let X m3iX (E) be the maximum depth of the shower produced by a 
primary of energy E. On average, the first interaction occurs at depth A, the mean free path 



of the initial particle. In this first interaction the initial particle splits into n(E) particles 
each carrying an average energy E/n(E). Therefore, we have 

X max (£) = A + X max (E/n(E)). (7) 

Assuming that X ma _ x (E) depends logarithmically on energy we get 

X max (£) = Alog 10 (E/n(E)) + B, (8) 

where A = XolnlO and B are constants. Xo = 37 gr/cm 2 is the electromagnetic radiation 
length. If we now assume that n(E) oc E A , we get 

X max (E) = A(l - A) log 10 (£) + B'. (9) 

This is the elongation rate theorem. We can now directly read the elongation rate from the 
above equation D\q = A(l — A). As stated previously, a change in the multiplicity growth 
with energy implies a change in the elongation rate. 

In Fig. ([I]) we show X max as a function of energy for the Fly's Eye and HIRES-MIA exper- 
iments. Data have been taken from references [|I],|J. The errors shown are only statistical. 
An additional systematic error of ~ 20 gr/cm 2 must be included in the data. The dash 
line represents our calculation for the slope parameter, D\o = 65 gr/cm 2 , (A = 0.24 from 
our simulations) for Fe-Air collision without fusion, normalized with the data at 6 x 10 17 
eV. The dotted curve has a slope parameter D w = 78 gr/cm 2 , which would imply a maxi- 
mum reduction in the slope of growing of multiplicities: from A ~ 0.24 to A ~ 0.09. The 
data from HIRES-MIA is not completely consistent with the Fly's Eye data. Statistical 
uncertainties are larger. The elongation rate obtained by the HIRES-MIA collaboration is 
very large, 95 gr/cm 2 . Notice that the elongation rate theorem predicts an elongation rate 
always less than 85 gr/cm 2 , an elongation rate larger would imply multiplicities decreasing 
with energy. Therefore, if the HIRES-MIA result is confirmed, a change of the composition 
is necessary to explain the data. In the figure we also show the energy region in which 
a phase transition is expected in Fe-Air collision using the string fusion model [19|. This 
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region corresponds to the range obtained from percolation theory, 1.1 < i] < 1.2, where rj is 
given by Eq.([3|). 

A word of caution is necessary in the use of the elongation rate theorem. The elongation 
rate theorem is based on the assumption that the energy is equally shared between the 
secondaries in the hadronic interaction. From this assumption immediately follows the 
logarithmic dependence of X max on the multiplicity. In realistics cases this assumption does 
not hold and one has to resort to simulations since no analytical formula is known for the 
X" max . We have parameterized for a number of models the dependence of X max on the 
change of multiplicity, see Ref. |19| for details. The results of a full Monte Carlo agree with 



our qualitative discussion. With this in mind we can conclude that to be able to explain 
the change on the elongation rate we need a change on the slope of the multiplicity from 
A ~ 0.24 to an essentially flat A ~ 0.09. The energy region for such change must be around 
5 x 10 17 eV. 

The experimental situation with the lateral distribution of muons p At (r) is clearer. Both 
in simulations and experiments it is found that the shape of the lateral distribution function 
for muons is rather independent of the primary's energy and composition. Therefore, at 
fixed distance to the core, r , Ph{tq) is proportional to the total number of produced muons 
in the shower. Under rather general arguments this number scales with energy |2l| 



P M (r ) oc = AE?, (10) 

where A is a normalization constant and /3 is the slope parameter. As mentioned previously, 
the slope parameter is found to be constant over a wide range of energies. This result is 
consistent with Yakutsk, Haverah Park pOj] , and with all the lower energy experiments. 

It is rather simple to calculate the slope parameter, (3, for a pionic cascade from a scaling 
argument similar to the elongation rate theorem. The number of muons is proportional to 
the number of charged pions at the maximum. The number of pions, at maximum, produced 
by a primary of energy E is given by 

N 7T (E ) = f n n [ dxP{x)N n {xE ), (11) 

J 
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where f n = 2/3 is the charged pion fraction, n the total pion multiplicity, and P(x) is the 
probability of producing a pion with a fraction of energy x of the primary energy. Assuming 
a scaling form, N n = AE@, we get 

l = Un f 1 dxP(x)x l3 = f n Z((3), (12) 
Jo 



where Z{(3) is the spectrum-weighted momentum For a given P(x) = 1/ndn/dx, the 
above equation gives an implicit equation for (3. It reduces to the textbook's expression if we 
assume energy equipartition, i.e. P(x) = S(x — l/n), which gives (3 = (l + log(/ 7r ))/log(n) ~ 
0.82, for n ~ 10. For realistic models the slope parameter (3 ranges between 0.7-0.9. In 
Eq.(|12"D the multiplicity enters explicitly in the left hand of the expression but also enters 
implicitly since the probability P(x) must verify total probability and energy conservation. 
Since Z((3) is a monotonically decreasing function of (3 for reasonable choices of dn/dx, a 
reduction of multiplicity induces a reduction on (3. Indeed this is what is observed for the 
dn/dx calculated for the model with and without fusion. The DPM model gives a value of 
(3 ~ 0.89 which agrees with that of the QGSJET model ||. In the presence of fusion the slope 
parameter is reduced and we get (3 ~ 0.72-0.77 depending on the specific implementation of 
the fusion model. This number was calculated both using Eq.QT2]) and by direct calculation 
with a Monte Carlo code. 

In Fig.(||), we show de density of muons at 600 m, p M (600), as a function of the energy for 
the AGASA measurements given as a parameterization and the HIRES-MIA measurements, 
shown as triangles. Also shown are the QGSJET results for pure iron and proton. Notice 
that the slope parameter measured by HIRES-MIA agrees with our slope parameter for the 
case of fusion. Our result, a change of slope parameter from that of pure iron to 0.77 at the 
energy where percolation is expected, is shown for comparison only. We can see that again 
it is consistent with the data. The slope measured by AGASA is different from the one 
calculated for either proton or iron for the QGSJET model. A rapid change on composition, 
as suggested by the HIRES-MIA data on A max , would imply a kink in the data for p Ai (600) 
at the same energy which is not seen. Instead, our result points towards a mild change on 
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the slope parameter, from 0.9-0.8 to 0.77 which would be hardly seen given the error bars 
in the data. 

There are a number of additional predictions in our scenario. The average pt in hadronic 
collisions should increase in the case of string fusion by about 10 - 20 %. This would produce 
flatter lateral distribution for the muon densities which could be observed. A particularly 
well-suited place to look for this effect would be in inclined showers. Inclined showers are 
composed esentially of high energy muons, and therefore are more sensitive to changes on 



the first hadronic interactions ||22|| . Given the current systematical and statistical errors we 
can not conclude that indeed cosmic ray experiments are observing a saturation of gluons 
or percolation of strings in hadronic interactions. High quality data with large statistics, 
as the expected from HIRES and the Pierre Auger observatories, are needed. However it is 
suggestive that all cosmic ray existing data are consistent with such interpretation. RHIC 
and LHC will measure the total multiplicity in the relevant energy region and ascertain 
whether a strong reduction of multiplicity takes place or not. But in any case, even if 
the change of composition is real, these effects must be taken into account in a complete 
simulation of cosmic ray showers. 
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FIGURES 




FIG. 1. Depth of shower maximum as a function of the logarithm of the primary energy as 
measured by Fly's Eye (circles) and HIRES-MIA (triangles). Full lines are fits to the data. Dashed 
and dotted lines are our prediction for a strong reduction on the multiplicity at an energy of 
~ 6 x 10 17 eV. Arrows mark the expected region where percolation occurs for Fe-Air collisions in 
the string fusion model. 
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FIG. 2. Logarithm of the muon density at 600 m as a function of the logarithm of the energy 
as measured by HIRES-MIA (triangles) and AG AS A (dashed line). Dotted lines are the errors 
of the AGASA parameterization. Also shown the prediction for the QGSJET model for pure iron 
(upper full line) and proton (lower full line) and our prediction for a change of slope as given in 
the text (marked). The arrows mark the position where percolation occurs for Fe-Air collisions in 
the string fusion model. 
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